Abstract
This notebook is a template for evaluating COVID-19 hospitalization forecast submissions from COVIDhub. After inputting a set of parameters (forecasters, COVID signals, etc), the template yields a comprehensive report on the predictions of COVID forecasters as well as their performance compared to the ground truth. The visualizations generated by the template offer an intuitive way to compare the accuracy of forecasters across all US states.
\[\\[.4in]\]
Every week, forecasters submit their hospitalization predictions to COVID-19 ForecastHub. In this report, we rely on an AWS bucket that contains the estimates of a handful of signals (e.g., COVID death, cases, hospitalization, etc). Furthermore, the AWS server stores an array of evaluation metrics of these forecasts (e.g., Absolute Error, Weighted Interval Score, and 80% Coverage). Alternatively, the data can be retrieved from the publicly accessible covidcast and covideval APIs.
s3a <- get_bucket("forecast-eval", region = "us-east-2")
s3b <- get_bucket("forecasting-team-data")
scores1 <- s3readRDS("score_cards_state_hospitalizations.rds", s3a,
region = "us-east-2")
scores1 <- subset(
scores1,
forecaster %in% union(params$forecasters, "COVIDhub-baseline"))
# The Hub always makes forecasts on Monday (dofw 2)
# We move all forecast_dates to the following Monday and shorten the ahead
wday_shift <- function(x, base_dofw = 2) (base_dofw - wday(x)) %% 7
scores1 <- scores1 %>%
mutate(ahead = ahead - wday_shift(forecast_date),
forecast_date = forecast_date + wday_shift(forecast_date))
# need to pass in the right forecaster here
scores <- s3readRDS(params$aws_scores, s3b) %>%
magrittr::extract2("reformatted.buggy.matched.scorecards") %>%
select(ahead, geo_value, forecaster, forecast_date, target_end_date,
actual, wis, ae, cov_80, value_50, data_source, signal,
incidence_period, direction_on_forecast_date)
direction <- scores %>%
select(forecast_date, direction_on_forecast_date) %>%
distinct()
# Logan's names are brutally long...
our_pred_dates <- unique(scores$forecast_date)
forecast_dates <- our_pred_dates
aheads <- unique(scores$ahead)
geo_values <- unique(scores$geo_value)
veggies <- c("asparagus", "broccoli", "chard", "daikon",
"escarole", "fennel", "garlic", "horseradish",
"jicama", "kohlrabi", "lettuce", "mushroom",
"nori", "okra", "parsnip", "radish", "squash",
"tomato", "watercress", "baseline1", "yuca")
names_table <- tibble(forecaster = unique(scores$forecaster), veggies = veggies)
names_table <- names_table[names_table$veggies %in% c("fennel", "jicama", "mushroom", "tomato", "watercress", "yuca"), ]
scores <- left_join(scores, names_table) %>%
filter(veggies %in% c("fennel", "jicama", "mushroom", "tomato", "watercress", "yuca")) %>% # small subset
select(-forecaster, -direction_on_forecast_date) %>%
rename(forecaster = veggies)
results <- scores1 %>%
select(ahead, geo_value, forecaster, forecast_date, target_end_date,
actual, wis, ae, cov_80, value_50, data_source, signal,
incidence_period) %>%
bind_rows(scores) %>%
filter(forecast_date %in% forecast_dates,
ahead %in% aheads,
geo_value %in% geo_values) %>%
left_join(direction) %>%
filter(!is.na(direction_on_forecast_date))
The target forecast dates are:
2020-11-16, 2020-11-23, 2020-11-30, 2020-12-07, 2020-12-14, 2020-12-21, 2020-12-28, 2021-01-04, 2021-01-11, 2021-01-18, 2021-01-25, 2021-02-01, 2021-02-15, 2021-02-22, 2021-03-01, 2021-03-08, 2021-03-15, 2021-03-22, 2021-03-29, 2021-04-05, 2021-04-12, 2021-04-19, 2021-05-03, 2021-05-10, 2021-05-17, 2021-05-24, 2021-05-31, 2021-06-07, 2021-06-14, 2021-06-21, 2021-06-28, 2021-07-05, 2021-07-12, 2021-07-19, 2021-07-26, 2021-08-02, 2021-08-09, 2021-08-16, 2021-08-23, 2021-08-30, 2021-09-06, 2021-09-13, 2021-09-20, 2021-09-27
The template will compile data of the following forecasters:
COVIDhub-ensemble.
For this analysis, all of Logan’s forecasters have been renamed:
kableExtra::kbl(names_table) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))
| forecaster | veggies |
|---|---|
| ______________________none oareigo zareino neimacpas | fennel |
| ______________________none zacougo dacouno neiwinnah | jicama |
| __________________never-up oareigo zareino neimacpas | mushroom |
| __nqcomb_logistic1_u/s/d/l oareigo zareino neimacpas | tomato |
| _marginal_cheating_u/s/d/l oareigo zareino neimacpas | watercress |
| marginal_logistic1_u/s/d/l oareigo zareino neimacpas | yuca |
\[\\[.07in]\]
Mean <- function(x) mean(x, na.rm = TRUE)
GeoMean <- function(x, offset = 0) exp(Mean(log(x + offset)))
facets.label = str_glue("{aheads} days ahead")
names(facets.label) = aheads
subtitle = sprintf("Forecasts made over %s to %s",
format(min(forecast_dates), "%B %d, %Y"),
format(max(forecast_dates), "%B %d, %Y"))
plot_wis_a <-
plot_canonical(results,
x = "ahead",
y = "wis",
aggr = GeoMean,
grp_vars = c("forecaster", "direction_on_forecast_date"),
facet_rows = "direction_on_forecast_date",
dots = TRUE,
base_forecaster = "COVIDhub-baseline") +
labs(title = subtitle,
x = "Days ahead",
y = "Geometric Mean WIS") +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
geom_hline(yintercept = 1, size = 1.5) +
scale_y_log10() +
scale_color_viridis_d() +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_wis_a, tooltip="text", height=800, width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))
plot_wis_a <-
plot_canonical(results,
x = "ahead",
y = "wis",
aggr = Mean,
grp_vars = c("forecaster", "direction_on_forecast_date"),
facet_rows = "direction_on_forecast_date",
dots = TRUE,
base_forecaster = "COVIDhub-baseline") +
labs(title = subtitle,
x = "Days ahead",
y = "Geometric Mean WIS") +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
geom_hline(yintercept = 1, size = 1.5) +
scale_y_log10() +
scale_color_viridis_d() +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_wis_a, tooltip="text", height=800, width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))
plot_cov80_a <-
plot_canonical(results,
x = "ahead",
y = "cov_80",
aggr = mean,
grp_vars = c("forecaster", "direction_on_forecast_date"),
facet_rows = "direction_on_forecast_date",
dots = TRUE) +
labs(title = subtitle, x= "Days ahead", y = "Mean Coverage 80") +
theme_bw() +
theme(plot.title = element_text(hjust = "center"),
legend.position = "bottom",
legend.title = element_blank()) +
scale_color_viridis_d() +
geom_hline(yintercept = 0.8, size = 1.5) +
guides(color = guide_legend(ncol = 2))
ggplotly(plot_cov80_a, tooltip="text", height=800, width=1000)
actuals <- results %>%
select(forecast_date, geo_value, actual) %>%
distinct()
traj_plot <- function(location) {
results %>%
filter(geo_value == location) %>%
ggplot() +
geom_line(data = actuals %>% filter(geo_value == location),
aes(forecast_date, actual)) +
geom_line(aes(target_end_date, value_50, group = forecast_date, color = forecaster)) +
facet_wrap(~forecaster, ncol = 2) +
theme_bw() +
theme(legend.position = "bottom") +
scale_color_brewer(palette = "Set1")
}
traj_plot("ny")
traj_plot("tx")
traj_plot("ca")
traj_plot("fl")
traj_plot("ga")
traj_plot("tn")
traj_plot("ms")
traj_plot("ut")
traj_plot("nh")